//本文件是gronber特殊高斯消去pthread并行化有关的各函数的合集
//注意，进行测试时可能需要分到多个文件中，由于函数来自多个工程，所以有一些变量的冲突什么的

/*
****我实现的算法有***
grobner					avx
grobner_pthread			avx_pthread
**********************
*/


#include<iostream>
#include<fstream>
#include<string>
#include<sstream>
#include<map>
#include<windows.h>
#include<tmmintrin.h>
#include<xmmintrin.h>
#include<emmintrin.h>
#include<pmmintrin.h>
#include<smmintrin.h>
#include<nmmintrin.h>
#include<immintrin.h>
#include<pthread.h>
#include<omp.h>
using namespace std;

#define NUM_THREADS 7

struct threadParam_t {
	int t_id;
	int num;
};

const int maxsize = 3000;
const int maxrow = 60000;
const int numBasis = 100000;

pthread_mutex_t lock;

long long head, tail, freq;

map<int, int*>iToBasis;
map<int, int*>ans;

int gRows[maxrow][maxsize];
int gBasis[numBasis][maxsize];


void grobner() {
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);

	int num = (flag == -1) ? maxrow : flag;
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	for (int i = 0; i < num; i++) {
		while (findfirst(i)!= -1) {
			int first =findfirst(i);
			if (iToBasis.find(first) != iToBasis.end()) {
				int* basis = iToBasis.find(first)->second;
				for (int j = 0; j < maxsize; j++) {
					gRows[i][j] = gRows[i][j] ^ basis[j];

				}
			}
			else {
				for (int j = 0; j < maxsize; j++) {
					gBasis[first][j] = gRows[i][j];
				}
				iToBasis.insert(pair<int, int*>(first, gBasis[first]));
				ans.insert(pair<int, int*>(first, gBasis[first]));
				break;
			}
		}
	}
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "Ordinary time:" << (tail - head) * 1000 / freq << "ms" << endl;
}

void avx_grobner() {
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);
	int num = (flag == -1) ? maxrow : flag;
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	for (int i = 0; i < num; i++) {
		while (findfirst(i) != -1) {
			int first = findfirst(i);
			if (iToBasis.find(first) != iToBasis.end()) {
				int* basis = iToBasis.find(first)->second;
				int j = 0;
				for (; j + 8 < maxsize; j += 8) {
					__m256i vij = _mm256_loadu_si256((__m256i*) & gRows[i][j]);
					__m256i vj = _mm256_loadu_si256((__m256i*) & basis[j]);
					__m256i vx = _mm256_xor_si256(vij, vj);
					_mm256_storeu_si256((__m256i*) & gRows[i][j], vx);
				}
				for (; j < maxsize; j++) {
					gRows[i][j] = gRows[i][j] ^ basis[j];
				}
			}
			else {
				int j = 0;
				for (; j + 8 < maxsize; j += 8) {
					__m256i vij = _mm256_loadu_si256((__m256i*) & gRows[i][j]);
					_mm256_storeu_si256((__m256i*) & gBasis[first][j], vij);
				}
				for (; j < maxsize; j++) {
					gBasis[first][j] = gRows[i][j];
				}
				iToBasis.insert(pair<int, int*>(first, gBasis[first]));
				ans.insert(pair<int, int*>(first, gBasis[first]));
				break;
			}
		}
	}
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "AVX time:" << (tail - head) * 1000 / freq << "ms" << endl;
}

void* grobner_pthreadfunc(void* param) {

	threadParam_t* p = (threadParam_t*)param;
	int t_id = p->t_id;
	int num = p->num;

	for (int i = t_id; i + NUM_THREADS < num; i += NUM_THREADS) {
		while (findfirst(i) != -1) {
			int first = findfirst(i);
			if (iToBasis.find(first) != iToBasis.end()) {
				int* basis = iToBasis.find(first)->second;
				for (int j = 0; j < maxsize; j++) {
					gRows[i][j] = gRows[i][j] ^ basis[j];

				}
			}
			else {
				pthread_mutex_lock(&lock);
				if (iToBasis.find(first) != iToBasis.end())
				{
					pthread_mutex_unlock(&lock);
					continue;
				}
				for (int j = 0; j < maxsize; j++) {
					gBasis[first][j] = gRows[i][j];
				}
				iToBasis.insert(pair<int, int*>(first, gBasis[first]));
				ans.insert(pair<int, int*>(first, gBasis[first]));
				pthread_mutex_unlock(&lock);
				break;
			}

		}
	}
	cout << t_id << "线程完毕" << endl;
	pthread_exit(NULL);
	return NULL;
}
void grobner_pthread() {
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);

	int num = (flag == -1) ? maxrow : flag;
	
	pthread_mutex_init(&lock, NULL);
	
	pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
	threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));

	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		param[t_id].t_id = t_id;
		param[t_id].num = num;
		pthread_create(&handle[t_id], NULL, grobner_pthreadfunc, &param[t_id]);
	}

	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		pthread_join(handle[t_id], NULL);
	}

	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "grobner_pthread time:" << (tail - head) * 1000 / freq << "ms" << endl;
	free(handle);
	free(param);
	pthread_mutex_destroy(&lock);
}

void* avx_pthreadfunc(void* param) {

	threadParam_t* p = (threadParam_t*)param;
	int t_id = p->t_id;
	int num = p->num;

	for (int i = t_id; i + NUM_THREADS < num; i += NUM_THREADS) {
		while (findfirst(i) != -1) {
			int first = findfirst(i);
			if (iToBasis.find(first) != iToBasis.end()) {
				int* basis = iToBasis.find(first)->second;
				int j = 0;
				for (; j + 8 < maxsize; j += 8) {
					__m256i vij = _mm256_loadu_si256((__m256i*) & gRows[i][j]);
					__m256i vj = _mm256_loadu_si256((__m256i*) & basis[j]);
					__m256i vx = _mm256_xor_si256(vij, vj);
					_mm256_storeu_si256((__m256i*) & gRows[i][j], vx);
				}
				for (; j < maxsize; j++) {
					gRows[i][j] = gRows[i][j] ^ basis[j];
				}
			}
			else {
				pthread_mutex_lock(&lock);
				if (iToBasis.find(first) != iToBasis.end())
				{
					pthread_mutex_unlock(&lock);
					continue;
				}
				int j = 0;
				for (; j + 8 < maxsize; j += 8) {
					__m256i vij = _mm256_loadu_si256((__m256i*) & gRows[i][j]);
					_mm256_storeu_si256((__m256i*) & gBasis[first][j], vij);
				}
				for (; j < maxsize; j++) {
					gBasis[first][j] = gRows[i][j];
				}
				iToBasis.insert(pair<int, int*>(first, gBasis[first]));
				ans.insert(pair<int, int*>(first, gBasis[first]));
				pthread_mutex_unlock(&lock);
				break;
			}
		}
	}
	cout << t_id << "线程完毕" << endl;
	pthread_exit(NULL);
	return NULL;
}
void avx_pthread() {
	int begin = 0;
	int flag;
	flag = readRowsFrom(begin);

	int num = (flag == -1) ? maxrow : flag;

	pthread_mutex_init(&lock, NULL);

	pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
	threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));

	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		param[t_id].t_id = t_id;
		param[t_id].num = num;
		pthread_create(&handle[t_id], NULL, avx_pthreadfunc, &param[t_id]);
	}

	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		pthread_join(handle[t_id], NULL);
	}

	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "avx_pthread time:" << (tail - head) * 1000 / freq << "ms" << endl;
	free(handle);
	free(param);
	pthread_mutex_destroy(&lock);
}



void reset() {
	memset(gRows, 0, sizeof(gRows));
	memset(gBasis, 0, sizeof(gBasis));
	RowFile.close();
	BasisFile.close();
	RowFile.open("被消元行.txt", ios::in | ios::out);
	BasisFile.open("消元子.txt", ios::in | ios::out);
	iToBasis.clear();
	ans.clear();
}

int readBasis() {          
	for (int i = 0; i < numBasis; i++) {
		if (BasisFile.eof()) {
			cout << "读取消元子" << i - 1 << "行" << endl;
			return i - 1;
		}
		string tmp;
		bool flag = false;
		int row = 0;
		getline(BasisFile, tmp);
		stringstream s(tmp);
		int pos;
		while (s >> pos) {
			if (!flag) {
				row = pos;
				flag = true;
				iToBasis.insert(pair<int, int*>(row, gBasis[row]));
			}
			int index = pos / 32;
			int offset = pos % 32;
			gBasis[row][index] = gBasis[row][index] | (1 << offset);
		}
		flag = false;
		row = 0;
	}
}

int readRowsFrom(int pos) {
	if (RowFile.is_open())
		RowFile.close();
	RowFile.open("被消元行.txt", ios::in | ios::out);
	memset(gRows, 0, sizeof(gRows));
	string line;
	for (int i = 0; i < pos; i++) {
		getline(RowFile, line);
	}
	for (int i = pos; i < pos + maxrow; i++) {
		int tmp;
		getline(RowFile, line);
		if (line.empty()) {
			cout << "读取被消元行 " << i << " 行" << endl;
			return i;
		}
		bool flag = false;
		stringstream s(line);
		while (s >> tmp) {
			int index = tmp / 32;
			int offset = tmp % 32;
			gRows[i - pos][index] = gRows[i - pos][index] | (1 << offset);
			flag = true;
		}
	}
	cout << "read max rows" << endl;
	return -1;
}

int findfirst(int row) {  
	int first;
	for (int i = maxsize - 1; i >= 0; i--) {
		if (gRows[row][i] == 0)
			continue;
		else {
			int pos = i * 32;
			int offset = 0;
			for (int k = 31; k >= 0; k--) {
				if (gRows[row][i] & (1 << k))
				{
					offset = k;
					break;
				}
			}
			first = pos + offset;
			return first;
		}
	}
	return -1;
}

void writeResult(ofstream& out) {
	for (auto it = ans.rbegin(); it != ans.rend(); it++) {
		int* result = it->second;
		int max = it->first / 32 + 1;
		for (int i = max; i >= 0; i--) {
			if (result[i] == 0)
				continue;
			int pos = i * 32;
			for (int k = 31; k >= 0; k--) {
				if (result[i] & (1 << k)) {
					out << k + pos << " ";
				}
			}
		}
		out << endl;
	}
}


